#------------------------------------------------------------------------------

#  读取消费过的会员信息：rawdata
data<-rawdata
data<-data[,sapply(data,class) %in% c('numeric','interger')]

# 数据初步整理

data<-data[outlier.delete(data),]

# 缺失值处理
data<-lossdelete(data)
lossrate(data)

# add.miss 自定义插补数值型变量
num
data<-add.miss(data,names(data)[lossrate(data)>0],num)

# 无法插补的缺失值去掉
data<-na.omit(data)

# 相关分析
library(psych)
mem.vcor<-corr.test(data,use="pairwise")
vnames<-names(data)[apply(mem.vcor,2,function(x) max(abs(i[i<1]))>0.4)]

n<-iif(length(vnames)>5,5,2)

data.svd<-svd(data[,vnames])
d<-data.svd$d[1:n]
u<-data.svd$u[,1:n]
v<-data.svd$v[,1:n]
data2<-u%*%diag(d)%*%t(v)

data<-cbind(data[-vnames],data2)

# 确定聚类个数
library(fpc)
pamk(data[sample(1:nrow(data),10000),],krange=5:15,criterion="multiasw",scaling=TRUE)$crit
pamk(data[sample(1:nrow(data),10000),],krange=5:15,criterion="ch",scaling=TRUE)$crit

pamk.result<-pamk(data[sample(1:nrow(data),10000),],krange=5:15,criterion="multiasw",scaling=TRUE)

clust<-slx.clust(scale(data),pamk.result)

dist(t(clust[[1]]$parameters$mean))








#==================================================================================
#==================================================================================
#   Appendix (Functions)
#==================================================================================
#==================================================================================

#---查看缺失率,删除缺失率高的字段
lossrate<-function(x){
  a<-apply(is.na(x),2,mean)
  a<-round(a,3)
  return(a)
}

lossdelete<-function(x){
  a<-apply(is.na(x),2,mean)
  a<-round(a,3)
  deletecols<-names(a[a>0.8])
  deleterows<-names(a[a>0&a<0.2])
  for(i in deletecols){
    x[,i]<-NULL
  }
  for(j in deleterows){
    x<-x[!is.na(x[,j]),]  
  }
  return(x)
}

#---插补缺失值
add.miss<-function(data,col,num){
  if(length(col)==1){
    i=col
    data[is.na(data[,i]),i]<-num
  }
  else{
    for(i in col){
      data[is.na(data[,i]),i]<-num[i]
    }
  }
  return(data)
}


# 训练模型

library(mclust)
slx.clust<-function(data,g){
  clust<-Mclust(data,G=g)
  res.anova<-apply(data,2,function(x) oneway.test(x~clust$classification)$p.value)
  res.mean<-apply(data[,which(res.anova<0.05)],2,function(x) tapply(x,clust$classification,mean,na.rm=T))
  return(list(clust,res.mean))
}

